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Acoustic wave fields propagating long ranges through 
the ocean are refracted by the inhomogeneities in the 
ocean's sound speed profile. Intuitively, for a given 
acoustic source frequency, the inhomogeneities become 
ineffective at refracting the field beyond a certain fine 
scale determined by the acoustic wavelength. On the 
other hand, ray methods are sensitive to infinitely fine 
features. Thus, it is possible to complicate arbitrarily the 
ray dynamics, and yet have the wave field propagate un- 
changed. This feature raises doubts about the ray/wave 
correspondence. Given the importance of various anal- 
yses relying on ray methods, a proper model should, at 
a minimum, exclude all of the fine structure that does 
not significantly alter the propagated wave field when 
the correspondence to the ray dynamics is integral. We 
develop a simple, efficient, smoothing technique to be 
applied to the inhomogeneities - a low pass filtering per- 
formed in the spatial domain - and give a characteri- 
zation of its necessary extent as a function of acoustic 
source frequency. We indicate how the smoothing im- 
proves the ray/wave correspondence, and show that the 
so-called "ray chaos" problem remains above a very low 
frequency (~ 15 — 25 Hz). 
PACS numbers: 43.30.Cq, 43.30.Ft, 43.20.Dk 

I. INTRODUCTION 

As acoustic waves propagate long ranges through the 
deep ocean, they are refracted by inhomogeneities in 
the ocean's sound speed profile. Roughly speaking, 
in the earth's mid-latitudes, temperature and pres- 
sure effectively combine to form a wave guide in the 
depth coordinate that vertically confines the prop- 
agation 1 . In addition to this overall structure, the 
ocean behaves as a weakly turbulent medium 2 that 
multiply scatters the acoustic waves mainly in the 
forward direction. Whether one is intrinsically in- 
terested in waves propagating through weak turbu- 
lence or in the state of the ocean through tomogra- 
phy 1 , ray methods are relied upon at various stages 
and levels of complexity in the resulting experimen- 
tal analyses 3-5 . It is therefore critical to understand 
the applicability and limits of these ray methods. 

Ray methods can only capture the physics of re- 
fraction and reflection, unless a geometric theory of 
diffraction is explicitly added 6 ' 7 . Intuitively, one ex- 
pects refractive effects to dominate diffraction when 



sound speed inhomogeneities are larger than the 
acoustic wavelength of the source. On the other 
hand, due to their pointlike nature, rays are sensitive 
to structures at all scales. Thus, one should be sus- 
picious of (non-diffractive) ray methods for models 
that have significant fine scale structure that are in- 
effective in refracting waves, but that fundamentally 
alter the rays themselves. Hence, certain fine scale 
structures in the model can be thought of as being 
physically irrelevant, i.e. having no influence on the 
wave propagation, and should be eliminated before 
applying a ray method analysis. The possibility of 
diffraction is very important, but should be dealt 
with separately and we do not discuss it further in 
this article. 

Another serious challenge for the applicability of 
ray methods that has been recognized in the past 
fifteen years or so, is the existence of ray chaos 8 ; see 
also earlier work in the field of quantum chaos . 
One typical argument goes that chaos introduces 
caustics, i.e. singularities, in ray methods at an ex- 
ponentially increasing rate with propagation time 
(range). Ray methods must therefore breakdown 
on a logarithmically short propagation scale, which 
renders them essentially useless. A significant body 
of work has shown that this need not be the case, 
and methods can be developed which are accurate 
to much longer propagation scales 13,14 . Even so, 
detailed ray methods tend to become rather bur- 
densome with the exponential proliferation of rays. 
Thus, resorting to statistical methods based on the 
chaotic properties of the rays is often attractive. 

These two reservations about ray methods, inclu- 
sion of physically irrelevant fine structures in the 
sound speed profile and ray chaos, have often been 
co-mingled. For example, it is possible to add very 
fine structure to a sound speed model that has no ef- 
fect on propagating waves and yet generates chaotic 
rays as unstable as one wishes. Our point of view 
is that the two issues should be disentangled, neces- 
sarily beginning with the removal of the physically 
irrelevant fine structures, whose characterization de- 
pends on the acoustic wavelength. We will come 
back to the ray chaos question, but leave a more de- 
tailed and complete analysis for follow-up work to 
this paper. 

Our purpose is, thus, to create a practical and 
easily implemented technique for smoothing inho- 
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mogeneities in a sound speed model, and to give 
prescriptions for the extent of smoothing needed as 
a function of source frequency. Toward these ends, 
it is not necessary to mimic a realistic ocean accu- 
rately with the model, but rather to include certain 
key features, such as a simple form for the waveguide 
confinement and the fluctuations due to the weak 
turbulence. It is more than sufficient to include scat- 
tering solely in the vertical spatial plane, to make 
the parabolic approximation 15 and to neglect larger 
mesoscale structures. A simple ocean sound speed 
model uses Munk's canonical model 16 to account for 
large scale effects due to temperature, pressure and 
salinity, and an efficient implementation scheme by 
Colosi and Brown 17 to generate much smaller inho- 
mogeneities due to the ocean's internal waves. Using 
their approach, the inhomogeneities have the statis- 
tics of the Garrett-Munk spectrum 18 . These fea- 
tures, though leading to a simplified model of the 
ocean, are more than adequate for investigating the 
length scale at which fluctuation features become 
important. Increased realism will be included in a 
future companion paper 19 . 

The outline is as follows. In Section II, the ocean 
sound speed model and the method for acoustic 
propagation are presented. The following section 
considers theoretical issues such as the convergence 
of the propagated wave field by asking the question: 
"does adding more modes in the internal wave ex- 
pansion cease altering the propagation beyond some 
maximum mode number?" . In Section IV, a smooth- 
ing of the expression for the internal wave sound 
speed model is introduced. This smoothing is very 
similar to the application of a low-pass filter - it re- 
moves most of the structures in the sound speed 
model that are shorter than a certain scale - but 
it is done directly in the spatial domain so that ray 
methods can easily be applied. Sensibly, the op- 
timal amount of smoothing necessary is a function 
of source frequency. We demonstrate the effects of 
smoothing on both the wave field propagation and 
on the phase space structures associated with the 
underlying ensemble of rays. This does not, in gen- 
eral, eliminate the consideration of ray chaos as the 
Lyapunov exponents are still positive, but it does 
remove a significant amount of the so-called "micro- 
folding" 20 of the phase space structures. We discuss 
how this can markedly improve the ray/ wave corre- 
spondence. 

II. THE ACOUSTIC PROPAGATION MODEL 

In a medium such as the ocean where density fluc- 
tuations are small, the wave equation accurately de- 
scribes the acoustic waves in which we are interested. 



The governing equation is 

^(r,t) = c 2 (r,t)VH(r,t) , (1) 

where Re{$(r, t)} is the acoustic pressure and c(r, t) 
is the sound speed at a location r and time t. The 
mean sound speed is roughly 1.5 km/s and if we 
consider a water parcel, the sound passes through 
it far faster than any variation in c(r, t) due to the 
internal waves; i.e. the internal waves travel several 
orders of magnitude more slowly than the acoustic 
waves. Hence, it is reasonable to 'freeze' the state of 
the ocean and let c(r , t) = c(r) . 

In anticipation of treating long range propagation, 
we assume that the scattering in the azimuthal direc- 
tion is negligible and the important components of 
the acoustic wave field propagation take place in two 
spatial dimensions with f = (z,r), where z is depth 
in the ocean and r is range from the source. Consider 
a constant frequency source, i.e a pure sinusoidal 
source of angular frequency u = 2nf with frequency 
/, whose amplitude is constant in time. Then, the 
wave field has a frequency response, $ u (z, r), where 
$(z,r,t) = ^ u (z,r) e~ iut . With this assumption, 
the wave equation reduces to the Helmholtz equa- 
tion in cylindrical coordinates 

V 2 $ w (z,r) + fc 2 (z,r)$ w (z,r) =0 , (2) 

where the wave number k(z, r) = u)/c(z, r). 

A. The Parabolic Equation 

For long range propagation, waves that propagate 
too steeply with respect to the horizontal strike the 
ocean bottom and are strongly attenuated. Since 
the surviving waves propagate at small angles with 
respect to the horizontal, a Fresnel approximation 15 
is possible which expresses the acoustic frequency 
response as the product of an outgoing cylindri- 
cal wave, e zk ° r /y/r and a slowly varying envelope 
function, * w (^,r), where the horizontal wavenum- 
ber ko w oj/c . Thus, 

gjfc (w)r 

*o,(«,r) = r) ^ . (3) 

Substituting Eq. (3) into Eq. (2) and dropping two 
small terms gives the parabolic equation 

^l:^^' 7 -) =-7^2^2 V I / -(^r) + U(z,r)^(z,r) . (4) 

Since the sound speed can be decomposed into 
the reference sound speed, c , and fluctuations, 5c, 
about the reference: c(z, r) = c + 5c(z, r) with 
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5c(z,r) « c , the potential is approximated as fol- 
lows: 



2 I \c(z,r) 




Sc(z, r) 



CO 



(5) 



In our calculations, we'll use the last form of Eq. (5) 
for simplicity. Notice that there is a direct analogy 
between this parabolic equation and the quantum 
mechanical Schrodinger equation through the sub- 
stitutions: t — > r, m — > 1, and h — > 1/fco- We use 
a symmetric split-operator, fast-Fourier-transform 
method to propagate the wave field 21 ' 22 ; see Ap- 
pendix A. 

The two terms neglected on the right side of 
Eq. (4) are 



8fcP^ (Z ' r)+ 2^^ (Z ' r) - 



(6) 



The first term falls off rapidly with range since 
k r >> 1. The second term is dropped due to the 
parabolic approximation which assumes that for a 
slowly varying envelope function, the curvature is 
weak. Note that up to this point, we have also 
dropped other terms from the propagation equation 
in assuming negligible azimuthal scattering and neg- 
ligible time dependence of the internal waves. See 
the discussion in Ref. 23 for more details on all of 
the terms that have been dropped and an order of 
magnitude estimate for the size of the various con- 
tributions. 



B. Ocean Sound Speed Model 

A simple model for the speed of sound in the ocean 
consists of two main components. The first compo- 
nent of the model is an adiabatic, large scale be- 
havior which is responsible for creating the ocean's 
'sound channel' - an effective wave guide for acous- 
tic propagation in the deep ocean. This general be- 
havior has a minimum sound speed at the sound 
channel axis, and varies slowly with latitude and 
season, with the sound channel axis moving toward 
the surface for higher latitudes and colder seasons. 
Mesoscale variability is neglected in this study. The 
second component of the model is local fluctuations 
in the sound speed due to the ocean's internal waves. 
These fluctuations are much smaller in magnitude 
than the wave guide confining behavior, but describe 
the range dependence. The model potential V(z,r) 
takes the form 

V( z r) = fe(z ' r ) = 5cw ^ z ) + ^iU^Z) ( ? ) 



where Sc wg represents the change of the sound speed 
due to the wave guide, which we take to be range in- 
dependent, and Sc iw represents the fluctuations due 
to internal waves. 



1. The Confinement /Wave Guide 

In the ocean, the main effects of pressure, temper- 
ature, and salinity create a minimum in the sound 
speed. Since sound bends toward regions of lower ve- 
locity, the shape of the sound speed profile refracts 
propagating waves toward the sound channel axis. 
This effect is captured in a smooth, average model 
proposed by Walter Munk 16 and is known as Munk's 
canonical model 
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where rj(z) — 2[z — z a ]/B, z a is the sound chan- 
nel axis, B is the thcrmoclinc depth scale giving the 
approximate width of the sound channel, and 7 is 
a constant representing the overall strength of the 
confinement. This model has its minimum speed at 
z — z a and captures the right exponential and linear 
trends near the surface and bottom. The parameters 
are chosen to be B = 1.0 km, z a = 1.0 km, c = 1.49 
km/s and 7 = 0.0113 km -1 , which are roughly con- 
sistent with the well known environmental measure- 
ments performed in the SLICE89 experiment 24 ' 25 . 



2. Internal Wave Sound Speed Fluctuations 

Internal wave fluctuations perturb the sound 
speed in the ocean through the resultant vertical 
motions of water parcels. They are responsible for 
multiple, weak, forward scattering of acoustic waves. 
A numerical scheme has been introduced by Colosi 
and Brown 17 , which allows efficient computation 
of a random ensemble of individual realizations of 
the typical sound speed fluctuations. This scheme 
conforms to the Garrett-Munk spectral and statis- 
tical phenomenological description of the internal 
waves 18 ' 26 and has the form 



~ = £ H e ^ fc " exp ) s MM(z)) , (9) 

3=1 k r V 



c 



CO 
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where we took £(z) = e~ z / B - e~ H / B with H the 
depth of the ocean. The prefactor e J - i fe r includes a 
random phase and magnitude factor for each j and 
k r in the sum; see Appendix B for further details and 
to infer a definition of & 3 \k r - Since the frequency of 
vertical motions lie between the inertial frequency, 
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due to the earth's rotation, and the buoyancy fre- 
quency, due to the local stratification, the sum over 
the horizontal wave vector k r has terms representing 
the superposition of internal waves with wavelengths 
in the range of 1 — 100 km. A maximum for the j- 
summation has been chosen as J m ax — 180, which 
gives structure down to the scale of roughly a meter. 
The modes, sin(j7r£(z)), are connected to the buoy- 
ancy profile which is assumed to have an exponential 
form. Although the form given in the text above for 
£(z) does not vanish precisely at the surface, it is 
sufficient for our purposes. 

The model should enforce that both the function 
Sc iw and its derivative vanish sufficiently smoothly 
at the surface. Thus, a surface filter is introduced 
which consists of multiplying Eq. (9) by a contin- 
uous function of depth with the properties that it 
vanishes above the ocean's surface, is unity below 
200 m, and has continuous first and second deriva- 
tives. In this way, 5ci W /co and its derivative vanish 
at the surface and are fully, smoothly restored be- 
low 200 m. Since the upper 200 m of the ocean can 
be quite complex with storms, seasonal fluctuations 
and latitudinal variability, there is no simple, general 
sound speed model near the surface; the surface fil- 
ter is adequate for our purposes. We will propagate 
waves for which very little energy will enter this re- 
gion, and thus, little effect of this surface smoothing 
will be relevant. The specific form we have chosen 
for the surface filter is 



g(z;z st ,T st ) 



for z' < -1/2 
h(z') for \z'\ < 1/2 

1 for z'> 1/2 



(10) 



where z' = (z — z st )/T st , the width is T st = 200 m, 
the center is z st — T s t/2 = 100 m, and the smooth 
function in between is 

h(z) = ^ + ^ S m(nz) + ^sm(3nz) . (11) 



C. Initial Wave Field 

The parabolic equation requires an initial wave 
field ^ u (z, r — 0) as input, which can then be propa- 
gated to the desired range of interest. It is important 
to understand the connection between the initial 
wave field and the localized, continuous wave source. 
Typical sources can be thought of as point sources 
whose acoustic energy disperses broadly. However, 
due to the previously mentioned fact that all the 
steeply propagating waves are strongly attenuated, 
we can instead propagate only that wave energy 
moving sufficiently close to the horizontal (within 
a spread of angles from the sound channel axis) that 



would have avoided the ocean's surface and bottom. 
Restricting the propagation angles limits the size of 
the vertical wave vector and necessarily creates "un- 
certainty" in the location of the point source. For 
our purposes, it is appropriate to choose ^ w (z,0) 
to be a minimum uncertainty wave packet. This im- 
plies using a normalized Gaussian wave packet of the 
form 



*u,(«,0) 



1 



2ttct2 



exp 



+ ik 0z (z - z ) 
(12) 



where zo centers the field, u z is the standard de- 
viation of the Gaussian intensity and kg z gives the 
propagating field an initial wavenumber in the z- 
direction. In all our calculations, we set k nz = 0, 
which maximizes the horizontally propagating en- 
ergy, and zq — z a , which centers the energy on the 
sound channel axis. 

A Fourier transform of Eq. (12) yields a complex 
Gaussian distribution of initial vertical wave num- 
bers, k z , centered at k 0z with standard deviation 
in intensity, at- Since a z and o\ are the variances 
of the intensity and not the amplitude of the wave, 
their relation is cr 2 — l/4a 2 . By a simple geometrical 
argument, a vertical wavenumber can be related to 
the horizontal wavenumber by k z — fc tan#, where 
9 is the angle with respect to the sound channel axis. 
In the next subsection, it is seen that p — tan# is 
a generalized momentum for a classical ray corre- 
sponding to the wave. Classical rays with the max- 
imum angle max just barely graze the surface or 
bottom, and thus, rays are limited in their verti- 
cal wave numbers. Yet, for Gaussian wave packets, 
all wave numbers are in principle present, though 
most are weighted negligibly by the tails. It is the 
width, <7fc, which determines if the wave contains 
wave numbers large enough for a substantial amount 
of the wave to hit the surface or the ocean floor. 
One can determine a proper Gaussian width, in or- 
der for only the Gaussian tails to pass the surface 
or bottom, in analogy with the limiting classical 
rays by letting the maximum classical wavenumber 
fco tan m ax correspond to three standard deviations 
out in the initial Gaussian wavenumber distribution, 
i.e. set 3<7/c = fcotan# max . Then 

(13) 



4/cq tan 2 m ax 



The explicit dependence of o z on the angular fre- 
quency, u), of the continuous wave source is realized 
using the approximate relation ko w uj/cq. 

The specific choice of 9 max depends on the verti- 
cal confinement. For the background confinement in 
Eq. (8), those rays departing the sound channel axis 
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with an angle of 8 = w/18 (10°) come within 150 m 
of the surface, and those with 9 = n/lb (12°) come 
within 30 m. The addition of internal waves to the 
sound speed model causes some rays to have a closer 
approach to the surface, so we will most often use 
Omax = 10° in this paper. 



D. The Classical Rays 

From the parabolic equation, one can derive a 
Hamiltonian system of equations for the position, 
z, and generalized momentum, p, of the collection 
of rays corresponding to the wave propagation. The 
Hamiltonian is given by H = p 2 /2 + V(z, r) and the 
potential is V(z, r) = Sc(z,r)/cQ. The equations are 



dz 
dr 
dp 
dr 



OH 

dp 
dH 

dz 



V 



dV(z,r) 
dz 



(14) 



Since dz/dr w Az/Ar = tan#, the generalized mo- 
mentum is p = tan#. The classical action T is cal- 
culated by imposing the initial condition To = and 
using the relationship 



dT _ dz 
dr ^ dr 



(15) 



Through the parabolic approximation, the classical 
action is directly related to the travel time, t, of the 
acoustic waves, where T — c r — r. 

The relevant rays to the wave propagation are 
those appropriate for a Gaussian wave packet 14,27 , 
which implies initial conditions in the neighborhood 
of (zo,Po). Since ko z — fcoPo = for the wave packet 
in Eq. (12), ray calculations are done in a neighbor- 
hood of po = 0. However, zo is taken to be on the 
sound channel axis, z a . 

The addition of range dependent internal wave ef- 
fects to the sound speed model causes the classical 
rays to be chaotic 8 . The stability matrix contains 
the information about whether the rays are unsta- 
ble (chaotic) or not 28 . At a fixed r, one has 



5p r 
Sz r 



) ( Sp a 

lr \ Sz 



where the stability matrix 
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qn qi2 

921 922 



/ dp. 




dp r 




zo 


dz 


dz r 


dz r 


\ dpo 


zo 


dz Q 



Po 



Pa 



Elements of this matrix evolve according to 
d 

——Cj r — K r Q r , 
dr 



(16) 



(17) 



(18) 



where Q r at r — is the identity matrix, and 



K r = 




(19) 



The system of differential equations Eqs. (14), (15) 
and (18), are solved using a 4th order Runge-Kutta 
method (where we have taken Ar = 100 m in all 
calculations) . 

The Lyapunov exponent, /i, is a measure of the 
rate at which the rays are deviating under small per- 
turbations. The relationship between the Lyapunov 
exponent and the matrix Q r comes through the trace 
(sum of the diagonal elements) of Q r , 



a= lim -ln|7Y(Q r )| . 



(20) 



If \Tr{Q r )\ grows exponentially, the Lyapunov ex- 
ponent is nonvanishing and positive, and the corre- 
sponding trajectory is chaotic. 



III. THEORETICAL CONSIDERATIONS 

Wave propagation should become increasingly in- 
sensitive to smooth perturbations as the scale of the 
perturbations decreases to the order of the smallest 
wavelength in the source and beyond. This issue 
does not arise in the horizontal coordinate of the in- 
ternal wave model in Eq. (9), since the fluctuation 
scales are much longer than the horizontal projec- 
tions of typical source wavelengths. However, this is 
an issue for the vertical fluctuations since Eq. (9) is a 
weighted superposition of a large number of vertical 
internal wave modes and presumably contains more 
detail than is necessary for accurate wave propaga- 
tion. There comes a point in the summation be- 
yond which the vertical modes begin to add physi- 
cally irrelevant features to the sound speed inhomo- 
geneities for a given source frequency. To determine 
the transition point where this occurs requires an 
understanding of the minimum wavelength structure 
in the propagating wave field, and an understanding 
of the power spectrum of individual vertical internal 
wave modes. The transition point, though, is not 
the only issue since higher modes contain a mix of 
physically relevant and irrelevant structures. These 
issues as well as their interplay are discussed here. 



A. The Vertical Mode Number Transition 

Intuitively, the vertical structures in the sound 
speed model responsible for refracting the wave are 
those that are larger than the minimum vertical 
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wavelength, \ m in, in the initial wave packet. Ex- 
pressions for X m i n can be obtained by using A = 
2ir/k and the previously noted geometrical relation 
k z = fco tan 9, 



Amm — 



2tt 



An 



c 



k Q tan 9 max tan 9 max f tan 9 



max 



(21) 



Recall that the ocean waveguide forces 9 ma x to be 
small so that the minimum vertical wavelength is 
always enhanced over the source wavelength, Ao- For 
Omax = 10°, this enhancement is roughly a factor of 
6. As a practical example, we note that some of the 
experiments conducted by the Acoustic Engineering 
Test (AET) 3 ' 4 use a broadband 75 Hz source. A 
pure 75 Hz source has a 20 m source wavelength. 
Thus, if the energy stripping due to the ocean surface 
and bottom is consistent with 9 max — 10°, then the 
wave propagation would have a minimum vertical 
wavelength scale of roughly 110 m. 

The vertical structures in the sound speed model 
in Eq. (9) arise through the superposition of vertical 
modes of the form e~ 3z / 2B sin(j7r(e~ z / B — e~ H l B )). 
Since the argument of the sine is nonlinear, each ver- 
tical mode contributes different oscillation lengths at 
different depths. The monotonicity of the argument 
illustrates that each mode has a "chirped" structure, 
i.e. each mode oscillates more and more rapidly as 
the surface is approached. To make this more pre- 
cise, an expansion of the argument of the sine re- 
veals that the local oscillation length as a function 
of depth is 



2Be z / B 



(22) 



Therefore, the j th internal wave mode contributes 
its shortest length contribution of 2B / j near the sur- 
face, with longer length scales at increasing depth. 
Each mode gives contributions to the sound speed 
fluctuations over a broad range of scales. 

Figure 1 illustrates the depth dependence and 
power spectrum of an internal wave mode. The 
power spectrum has a fairly sharp high frequency 
(short length scale) cutoff from the structures added 
near the ocean surface and a slowly decaying tail 
for the lower frequencies (longer length scales). The 
broad tail for an individual mode indicates that 
many different modes contribute to a particular size 
feature in the internal wave model. 

The transition vertical mode number 3 trans can be 
identified as that point where the vertical modes be- 
gin to introduce structure smaller than \ m in- Thus, 
setting Eqs. (21) and (22) equal to each other and 
solving for j gives 



Jtrans — 
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2Bt&n9 max _ 2Bft&n9 max 



Ao 



co 



(23) 
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FIG. 1. Example of the single vertical internal wave 
mode for j — 25. The upper plot illustrates its depth 
dependence, 

Vj = g(z;z st ,T st ) e-' Az ' 2B sm(jn(e-^ B - e - H ^ B )), where 
g is the surface filter defined in Eq. (10) and the lower 
plot is the power spectrum, P, of Vj. 



The calculation of Jtrans does not reflect that each 
vertical mode is weighted in Eq. (9) by the coeffi- 
cients e J - i fe r , which we numerically found to have root 

mean square decay \J^2k r \ e j-k r \ 2 ~ j 1 ' 1 f° r large 
j. Thus, the higher vertical modes have a slowly 
decreasing weighting. The acid test of the effects of 
both the diminishing amplitudes and the detectabil- 
ity of features by the wave is to look at the sensitiv- 
ity of the wave field to variations in the value for the 
j-summation cutoff in Eq. (9). 

B. Wave Field Convergence 

We can investigate the convergence of the wave 
field propagation by using different values for the 
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j-summation cutoff in Eq. (9) to generate various 
sound speed media. The value of the cutoff leading 
to a converged wave field, denoted by J w , is the min- 
imum number such that by including higher modes 
there is relatively little change in the wave propaga- 
tion. We do not have a simple intuitive argument 
that gives an expression for J w , but instead rely on 
numerical simulations to determine reasonable val- 
ues. 

In order to discuss quantitatively what is meant 
by 'little change to the wave propagation', it is nec- 
essary to have a measure of the similarity of two 
wave fields. An ideal measure is the absolute value 
squared of the overlap (inner product). For two 
sound speed potentials that differ by AV, the quan- 
tity Cav is defined as 



C A v(r) 



dz ^ Ay (z,r)* w (z,r) 



(24) 



where ^ u (z,r) is understood to be the wave field 
propagated to range r with the full potential and 
\t* AV (z, r) is the same initial state propagated using 
the potential which differs from the full potential by 
AV. It is convenient to normalize the propagating 
wave fields to unity since this is preserved under the 
unitary propagation of the parabolic equation. With 
this choice, the measure gives unity only if the two 
wave fields are identical. The greater the reduction 
from unity, the greater the differences between the 
two propagations, i.e. the lower the faithfulness or 
fidelity of the propagations. 

The full wave propagation is compared to wave 
propagation for various values of mode number cut- 
off J < J max- Thus, AV is the internal wave sum 
for j in the interval [ J + 1 , J ma J ■ Since deviations 
of C\v{r) from unity indicate that features in the 
modes [J + l,J ma x] were detectable by the wave 
propagation, the value of J where CAv{ r ) breaks 
appreciably from unity determines J w . 

Sound waves with source frequencies of 25, 75, 150, 
and 250 Hz were propagated to r = 1000 km; sec 
Appendix A for details regarding the propagation. 
Figure 2 demonstrates the dependence of Cav i r ) on 
J. To interpret this figure, consider the curve for 75 
Hz. Above J = 50, CAv( r ) > -99 and remains near 
unity. Thus, we can say that here J w w 50. Using 
higher internal wave modes adds no more realism, 
and only slows down the calculations. A similar ar- 
gument for the other frequencies gives the values of 
J w listed in Table I. Note that J w increases just a lit- 
tle more slowly than linear in source frequency due, 
in part, to the decreasing weightings. 




FIG. 2. G'Av(r) as a function of J for the source fre- 
quencies of 25, 75, 150 and 250 Hz (corresponding to the 
curves from left to right, respectively) at a range of 1000 
km. 



TABLE I. Comparison of key parameters for a few 
viable long range propagation frequencies. Both J„ and 
\° pt (see the next section) were determined using a con- 
servative 0.99 criterion for the value of the Cav at 1000 
km in Figs. 2 and 6. Other choices for the criterion, 
propagation range, etc... could lead to somewhat greater 
differences; however the dependences are rather weak. 
For each calculation, 6 ma x = 10°. Note the minimum 
wavelength feature, A min , in the initial wave packet is 
extremely close to X° pt . 



Frequency 


Jtrans 




Ao 


Xmin 




(Hz) " 






(km) 


(km) 


(km) 


25 


6 


20 


0.060 


0.340 


0.308 


75 


18 


50 


0.020 


0.113 


0.106 


150 


36 


90 


0.010 


0.056 


0.060 


250 


60 


145 


0.006 


0.034 


0.034 
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Since Cav is inherently range dependent, deter- 
mining J u from a plot of Cav for a single range is po- 
tentially inappropriate for longer ranges. Yet, Fig- 
ure 3 illustrates that the range dependence of Cav 
is fairly weak for a frequency of 75 Hz. Increasing 
the range from 1000 km to 4000 km for J u = 50 
decreases Cav very little from 0.99 to 0.96. Thus, 
J w = 50 is a conservative choice even for ranges up 
to 4000 km. 




FIG. 3. Cav(t) as a function of J for the ranges 
of 1000, 2000, 3000 and 4000 km (corresponding to the 
curves from left to right, respectively) for a frequency of 
75 Hz. 

For reasonable source frequencies, it is clear that 
even an optimal choice for J w will leave a significant 
amount of oscillations in the model on a scale much 
smaller than A mi „ (since J u » Jtrans)- Presum- 
ably, these oscillations are physically irrelevant for 
the wave propagation, but it is worthwhile study- 
ing more precisely where the cutoff lies within the 
context of long range propagation. 

IV. FILTERING THE PHYSICALLY 
IRRELEVANT FEATURES 

Since we have taken a smooth background sound 
speed model, the physically irrelevant features of the 
sound speed model can be removed by filtering the 
high frequency components from the internal wave 
sound speed model, 5ci w (z,r). The ideal approach 
would be through the application of a low pass fil- 
ter: Fourier transform 6ci w (z, r) for a fixed range to 
a frequency domain, apply a filter that removes the 
high frequencies and Fourier transform back to give 
the physically relevant portion of 5ci w (z,r). There 
are several drawbacks with respect to proceeding 
this way. The Fourier transforming back and forth 
is computationally expensive, creates a problematic 



ocean surface, and severely complicates the ray cor- 
respondence; the same would be true using a convo- 
lution technique. Instead, we develop a smoothing 
that can be directly applied to each vertical mode in 
the spatial z domain and serves as a very good ap- 
proximation to a low-pass filtering in the frequency 
domain. It takes advantage of the monotonicity of 
the chirped structure of the individual internal wave 
modes. The spatial filtering method simplifies the 
ray equations enormously and allows first and sec- 
ond derivatives to be evaluated exactly, as opposed 
to numerically, which is an unstable operation. 

A. The Smoothing 

Due to the precise oscillatory nature of each verti- 
cal mode, a good approximation to a low-pass filter 
can be accomplished by removing the upper portion 
of the vertical mode that contains oscillations that 
are smaller than the smoothing parameter, A s . This 
involves multiplying each vertical mode by the func- 
tion g(z; z sm ,T sm ) defined in Eq. (10). This filter 
is centered at the depth such that the local length 
scale is A s , which by inversion of Eq. (22), gives the 
mode-dependent depth z sm = Bln(j\ s /2B). Note 
that j must exceed 2B/X S in order for the filter to be 
below the ocean surface, which is where it begins to 
have an effect. This is consistent with the shortest 
length contribution of each mode being 2B j j at the 
surface. The width of the filter is carefully chosen to 
be T sm — 2. 0A S so that it does not cut off too sharply 
thereby introducing high frequency components into 
the model. If the width were chosen much greater, 
amplitudes of physically relevant length scales would 
be reduced. 




0.1 0.2 0.3 0.4 



X (km) 

FIG. 4. Effects of the amount of smoothing on the 
power spectrum, P, of Sc iw /co- The dashed line is the 
power spectrum of the unsmoothed full potential and 
the solid line is the power spectrum of the smoothed full 
potential for A s = 0.20 km. 
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Figure 4 shows the power spectrum of the sound 
speed model with and without smoothing; it is il- 
lustrated with a value, A s = 0.2 km. The power 
spectrum remains relatively unchanged for length 
scales greater than 0.2 km, but the length scales 
below 0.2 km are significantly dampened out of the 
model. This is evidence that a smoothing parameter 
of A s = 0.2 km is doing exactly what it was designed 
to do: it is filtering out features on scales below 0.2 
km, but leaving features above 0.2 km in the model. 
Figure 5 shows the smoothed sound speed potential 
and the portion of the potential, AV, filtered by the 
smoothing. It is clear from these figures that the os- 
cillations in the unsmoothed potential which have a 
length scale of less than 0.2 km have been removed, 
while larger oscillations have been preserved. 



B. Estimating the Optimal Smoothing 
Parameter 

The optimal smoothing parameter, \° pt , would be 
such that only those features in the model that are 
not detectable by the wave would be removed. Intu- 
itively, A° pt would be very close to A mi „ of Eq. (21). 
In order to test this intuition, we again use Cav{^) 
defined in Eq. (24), where here AV is the high fre- 
quency portion of the internal wave sum, which the 
smoothing removes, and the other potential is the 
full unsmoothed sound speed model. X° pt is deter- 
mined to be the maximum value of A s up to which 
Cav remains nearly unity but deviates significantly 
beyond. 
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FIG. 5. Effects of the amount of smoothing on the 
full potential, Sc im /co- In the upper panel, the dashed 
line is the unsmoothed potential and the solid line is 
the smoothed potential for A s = 0.20 km. In the lower 
panel, the difference, AV, between the smoothed and 
unsmoothed potential is displayed. 
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FIG. 6. Cav(t) as a function of A s for source fre- 
quencies of 25, 75, 150 and 250 Hz (corresponding to the 
curves from right to left respectively) at a range of 1000 
km. 



As in the previous section, source frequencies of 
25, 75, 150, and 250 Hz were propagated to r = 1000 
km with J chosen for each frequency to be that value 
of J w in Table I. Figure 6 demonstrates the de- 
pendence of CAv(r) on different values of A s and 
its interpretation is similar to that done for Fig. 2. 
Consider the curve for 75 Hz. Above A s w 0.1 
km, Cav breaks significantly from unity giving the 
optimal smoothing of the sound speed model for 
a 75 Hz source to be X opt w 0.1 km. Smooth- 
ing less than this allows high frequency features to 
remain in the model which have no effect on the 
wave propagation. Table I summarizes the results 
which all agree closely with the intuitive idea that 
A° P * ~ X m in = Ao/ tan# mox . 

For a fixed A s , the higher source frequencies lead 
to a reduced value of Cav- This indicates that 
the high frequency components of AV are more de- 
tectable by a high frequency source than by a low 
frequency source. This fully supports the age-old in- 
tuitive concept that high frequency waves can detect 
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smaller features than low frequency waves, and that 
the appropriate detection scale is a wavelength. A 
long range propagation experiment utilizing a source 
frequency / only detects that portion of the internal 
wave power spectrum with features longer than the 
length scale \ min = c Q / ftan6 max . 

C. Effects of Smoothing on Phase Space 
Structures 

Classical ray methods can be used to reconstruct 
propagating wave fields in detail through the use of 
semiclassical Green functions 29 . The semiclassical 
approximation to the wave field is 

* sc (z,r;fc ) = *^2Aj(z,r)exp[ikoTj(z,r) - invj/2] , 

3 

(25) 

where the sum is over all ray paths labeled by j that 
begin at the source and end at a depth z for a given 
range r. The phase contribution of a path is related 
to its classical action, Tj, the source wavenumber, 
fco, and the number of caustics, Vj. The amplitude 
contribution of a path, Aj, is related to its stability 
matrix elements; see Ref. 30 for a readable account. 
This discrete set of paths becomes continuous if wc 
consider all z. Thus, there is a continuous set of rays 
that underlies the full construction of ^f 8C (z,r;ko) 
at a given range. A powerful analysis of the prop- 
erties of this set comes by considering the rays in 
the phase space formed by all allowable points given 
by position and conjugate momentum. Viewed in 
phase space, the continuous set of rays underlying 
the wave field propagation (in the single degree of 
freedom problem being discussed here) forms a con- 
tinuous, self-avoiding line which is called a manifold. 
As the range increases, the manifold evolves into a 
rather wild-looking "spaghetti" . The more chaotic 
the system, the wilder the appearance of the mani- 
fold. 

The construction of Eq. (25) relics on the use of 
stationary phase approximations, which can only be 
applied reliably when the phase between successive 
stationary phase points is greater than order unity. 
Care must be taken in defining the meaning of suc- 
cessive in this context. Thus, Eq. (25) breaks down 
when [T 3 (z,r) - T jr (z,r)] < k^ 1 = \ /2ir where j 
and j' are the classical paths/rays corresponding to 
successive stationary phase points. We term this the 
'area-(Ao/27r) rule' (the translation to this problem 
of the area-fi. rule of Refs. 9 and 10). See Refs. 13 
and 14 for a detailed presentation of the breakdown 
of the stationary phase approximation in quantum 
chaotic systems. 



The breakdown of stationary phase is intimately 
related to how the manifold winds and folds its way 
through phase space. The difference in the classical 
action for two rays is related to the areas in phase 
space between the folds of the evolving manifold and 
the vertical line of the final depth, z, whose inter- 
sections with the manifold specify the rays. If these 
areas become smaller than Xq/2tt, then stationary 
phase breaks down for that pair of rays and we say 
that the two stationary phase points are coalescing. 
By drawing the manifold and filling in areas of Xo/2n 
in the folds, one can immediately see where prob- 
lems, such as caustics which produce infinite ampli- 
tudes, will be occurring in the semiclassical construc- 
tion. In the simplest case of two coalescing points, an 
Airy function uniformization is possible. However, if 
so many coalescing pairs occur that they cannot be 
isolated from each other, uniformization effectively 
is no longer possible, and the semiclassical approxi- 
mation has broken down. 

In the work of Simmen, Flatte, and Wang 20 , they 
show how the fine features in the internal wave field 
lead to a phenomenon they termed "micro-folding" 
in which tiny folds arc densely found along the man- 
ifold. Clearly, for typical source frequencies in long 
range propagation, the neighborhoods of the micro- 
folds violating the area-(Ao/27r) rule overlap every- 
where with each other. Thus, one anticipates a dense 
set of singularities in the semiclassical approxima- 
tion rendering the approach useless. 
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FIG. 7. Smoothed phase space manifold. The solid 
line is the phase space plot for a dense set of trajectories 
with launch angle € [4°,8°] propagated for 50 km in 
the unsmoothed ocean model. The dashed line is for the 
same set of trajectories, but for a smoothing parameter 
of A s = 0.10 km. All the trajectories started on the 
sound channel axis. The hatched rectangle is a reference 
area for physically irrelevant microfolds and has an area 
Xo/2-k, which corresponds to a 75 Hz source. 

Herein lies the advantage of smoothing the ocean 
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sound speed model ol physically irrelevant features 
before making the ray correspondence. Presum- 
ably the bulk of the micro-folding is related to 
fine features which are ignored by the wave prop- 
agation. The smoothed system contains only that 
structure necessary to describe the wave propaga- 
tion so it should have fewer micro-folds. Figure 7 
illustrates the effects of smoothing on a set of tra- 
jectories. One can see that the smoothed mani- 
fold tracks the unsmoothed manifold along its length 
very well. A more detailed example of micro-folding 
is illustrated in Fig. 8 for a range of 100 km. Notice 
how the smoothed manifold completely eliminates 
this particular micro-folded structure for a smooth- 
ing parameter of A s = 0.1 km (appropriate for 75 
Hz). Eleven, non- isolated pairs of coalescing station- 
ary phase points were eliminated by the smoothing. 
Only a well behaved piece of the manifold with no 
coalescing pairs remains. Thus, there are fewer loca- 
tions leading to singularities and breakdown in the 
semiclassical construction for the smoothed system, 
yet it is describing the same propagated wave. We 
leave the full semiclassical reconstruction for future 
work. 
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FIG. 8. Smoothing of the micro-folds. The solid line 
is the phase space plot for a dense set of trajectories 
with launch angle 9 G [7°,8°] propagated for 100 km in 
the unsmoothed ocean model. The dashed line is for the 
same set of trajectories, but for a smoothing parameter 
of A 3 = 0.10 km. All the trajectories started on the 
sound channel axis. The hatched rectangle is a reference 
area for physically irrelevant microfolds and has an area 
Xo/2-k, corresponding to a 75 Hz source. 



D. Effects of Smoothing on Lyapunov Exponent 

The following question naturally poses itself from 
the results of the previous section, "if smoothing the 
inhomogeneities reduces the number of folds, per- 
haps it is eliminating the ray chaos that was dis- 



covered in Ref. 8?" This turns out not to be the 
case. The Lyapunov exponents for smoothed sys- 
tems do not vanish. The Lyapunov exponent, fj,, as 
defined in Eq. (20), requires the infinite range limit, 
which due to the maximum range of the ocean, is 
not very sensible. Instead, it is much more relevant 
to work with finite-range Lyapunov exponents 31,32 . 
The stability matrix, Q r , as defined in Eq. (17), is 
calculated for a classical ray starting on the sound 
channel axis with an initial angle 9 and propagated 
for a range r. If \Tr Q r \ is growing exponentially 
with range, then the ray is unstable or chaotic and 
the following relationship can be inverted to obtain 
the finite-range Lyapunov exponent 

\Tr Q r \ = e» r + e-» r . (26) 

Excluding a few highly abstract systems, this /j, fluc- 
tuates as a function of range and from one ray to 
the next. In fact, for typical chaotic systems and 
the internal wave problem here, \Tr Q r \ is close to 
being lognormally distributed, or from a different 
point of view, the finite-range Lyapunov exponents 
give something close to a Gaussian density 31,32 . The 
finite-range Lyapunov exponents are launch angle 
dependent 33 . Figure 9 shows histograms of the 
finite-range Lyapunov exponents for a range of 1000 
km for a range of ray angles. 
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FIG. 9. Probability distribution of finite-range Lya- 
punov exponents. The range of propagation is 1000 km 
and each probability distribution consists of 4,000 tra- 
jectories within a uniform distribution of launch angles. 
For the solid line, |0| £ [0°,2°], for the dashed line, 
|0| € [4°,6°], and for the dotted line, |0| G [8°, 10°]. 
Each probability distribution was obtained by averaging 
over a Gaussian window of the corresponding histogram. 
The smoothing parameter is X s = 0.10 km and all the 
trajectories started on the sound channel axis. 

It turns out that the mean of the finite-range Lya- 
punov exponents is the usual infinite-limit Lyapunov 
exponent (as long as one has propagated beyond a 
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transient range of a few Lyapunov lengths). Letting 
the brackets <> denote averaging over many rays, 



A*o = - < ln|7Y Q r \ > ■ 
r 



(27) 



If one averages before taking the natural logarithm, 
one gets a second stability exponent which is not the 
Lyapunov exponent, but rather a related one: 



fj, 



1 

27 



ln(< \Tr Q r \ 2 >) ■ 



(28) 



The relationship between /iq and \x for a Gaussian 
density is through the variance of the distribution of 
the finite-range Lyapunov exponents 



2 M - Mo 



(29) 



These two stability exponents fix the Gaussian den- 
sity completely. Figure 10 illustrates the dependence 
of yu > fi an d the distribution on the smoothing pa- 
rameter A s . 
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FIG. 10. Average Lyapunov exponents, fxo and p,, and 
probability distribution as a function of the smoothing 
parameter, A s . Both plots are for a range of propaga- 
tion of 1000 km. The upper plot is an average of 2,000 
trajectories within a uniform distribution of launch an- 
gle 9 e [—10°, 10°]. The solid line is no and the dashed 
line is p,. The lower plot is the same as the previous 
plot except that \9\ € [8°, 10°] in both curves and the 
smoothing parameter is varied. The solid line is for a 
smoothing parameter of A s = 0.10 km and the dashed 
line is for a smoothing parameter of \ s — 0.30 km. A 
narrow peak near the origin exists in the dashed curve, 
which indicates a non-negligible fraction of stable trajec- 
tories. 



Although, there is still ray chaos, the Lyapunov 
exponent is monotonically decreasing with increased 
smoothing, but unless a smoothing greater than 0.10 
km is applied, /iq does not decrease appreciably. At 
some point, beyond a smoothing somewhere in the 
neighborhood of 0.3 — 0.5 km, a large fraction of 
the rays behave stably. Note in Fig. 10 that for 
A s = 0.30 km, a significant fraction of the rays have 
become stable, i.e., they have a Lyapunov exponent 
equal to zero. Using the relation between frequency 
and optimal smoothing, for source frequencies in the 
neighborhood of 15 — 25 Hz, there is a transition be- 
low which the ray chaos problem due to the inter- 
nal wave inhomogeneities effectively disappears and 
above which it remains important over ocean basin 
scale propagation ranges. Though the background 
profile used for this study is somewhat simplistic, 
surprisingly these results seem to be consistent with 
some very low frequency experiments. In particular, 
data from the Alternate Source Test (AST) clearly 
shows that 28 Hz receptions have a more stable ar- 
rival pattern than the 84 Hz receptions for transmis- 
sion over a 5000 km range 34 . 



V. DISCUSSION 

In probing the state of the ocean, it is impor- 
tant to understand what information is carried in 
the wave propagation for a given source frequency. 
Intuitively, fluctuations in the ocean sound speed on 
scales shorter than an acoustic wavelength should 
be ineffective sources of refraction for a sound wave 
in the ocean. Though, parabolic equation simula- 
tions arc unaffected by the inclusion of physically 
irrelevant fine scale fluctuations in the sound speed 
model (except for the resulting slower computation 
time), this inclusion worsens the correspondence of 
ray methods to the wave propagation. On the other 
hand, ray methods are sensitive to infinitely fine 
scale structures. Those fine structures that are not 
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detectable by the wave propagation lead to physi- 
cally irrelevant micro-folds in the phase space man- 
ifolds for the rays. These folds lead to unwanted 
singularities and the breakdown of semiclassical ap- 
proximations. Smoothing of the internal wave sound 
speed model allows a significant reduction in the ex- 
tent of micro-folding and this must lead to a better 
ray/wave correspondence. 

In our study, we noted that the chirped struc- 
ture of each of the internal wave modes gives con- 
tributions to the sound speed fluctuations over a 
broad range of scales. Thus, limiting the number 
of vertical modes used in an internal wave sound 
speed model does not fully resolve the issue of phys- 
ically irrelevant fine structure. For the specific con- 
struction of Colosi and Brown, our calculations gave 
frequency dependent values for the number of ver- 
tical modes J w necessary in the model. For fre- 
quencies of {25, 75, 150, 250} Hz, we found that the 
wave field propagation is essentially converged for 
J u = {20, 50, 90, 145}, respectively. However, for the 
same set of frequencies, modes greater than the tran- 
sition modes Jtrans = {6,18,36,60}, respectively, 
add structures on a finer scale than A m j n . Hence, 
each mode contains a large spread of frequency con- 
tributions so that a low-pass filtering of each vertical 
mode is needed. 

In order to remove physically irrelevant structures, 
we constructed an approximation in the position do- 
main to a low-pass filter by taking advantage of 
the monotonicity of the chirped structure of each 
mode. The accuracy of this approximation (though 
not shown in this paper) was very good for indi- 
vidual modes. The spatial filtering method that we 
developed gives three main advantages: reducing re- 
quired computations, better behavior in the neigh- 
borhood of the ocean's surface, and simplicity with 
respect to making the ray correspondence. With this 
study, it was found that the vertical scale of interest 
for the vertical fluctuations is not the source wave- 
length, Aq = Co//, but rather the minimum verti- 
cal wavelength present in the wave field, which con- 
tains the additional projection factor (tan# TOax ) ; 
see Eq. (21). 9 max is the largest angle with respect 
to the horizontal that waves can propagate with- 
out being stripped out by bottom interactions and 
is typically in the neighborhood of 10° — 12° in the 
ocean's mid-latitudes. For these values of 6 max the 
minimum vertical wavelength is roughly 5 — 6 times 
Ao; i.e. relevant vertical structures are much larger 
than that implied by Aq. 

Additionally, from the results in Table I, J w scales 
more slowly with increasing frequency than Jtrans- 
This appears to be due to the decreasing weighting 
of the terms in Eq. (9) , which directly influences the 
convergence of the wave field propagation and the 



value of J u . If this trend were to continue, then at 
a sufficiently high frequency, Jtrans would overtake 
J u in value. Beyond this frequency, low-pass filter- 
ing would no longer serve any purpose; one could 
simply choose an appropriate J u . We do not at- 
tempt to extrapolate to that point here using our 
calculations and model, but note that wherever it is, 
the frequency would be so high that very long-range 
acoustic propagation would not be possible due to 
dissipation. However, it may be useful in the con- 
text of short range acoustic experiments using much 
higher frequencies to establish a cross-over frequency 
with a more realistic model. 

We found that smoothing the internal wave sound 
speed fluctuations docs not, in general, eliminate 
the problems associated with ray chaos. The Lya- 
punov exponents are positive and significant unless 
the smoothing scale exceeds 300 - 500 m. Thus, in 
this simplified model, ray chaos continues to be an 
important issue for source frequencies above the 15 
- 25 Hz range. 

A number of difficulties arise in the study of 
chaotic systems. For example, the exponential pro- 
lification of rays, makes it impractical to carry out 
ray methods. A common technique to overcome 
these difficulties is to apply various statistical meth- 
ods whose justification derives from the chaos itself. 
However, even if you wish to apply these statistical 
methods, the validity of semiclassics is still an issue. 

Though it is known in the literature that the back- 
ground sound speed profile can dramatically affect 
the complexity of the ray dynamics, it is still a ques- 
tion for investigation as to how significant these ef- 
fects are on the wave propagation. Here we use 
Munk's canonical model as a simple, smooth back- 
ground profile, which is sufficient for a study of the 
removal of physically irrelevant structures. However, 
before inferring detailed properties of long range ex- 
perimental data, it would be good to employ a more 
realistic background sound speed wave guide. In 
fact, this would require a method for removing fine 
scale structures from the background in addition to 
the internal wave model and would not likely be sub- 
ject to as simple a spatial filtering scheme as we used 
for the internal waves. We will address these issue 
in a forthcoming paper. 

A number of interesting questions remain or 
emerge from our results. Our computations did not 
use pulsed sources, which can be expressed as an 
integral over a range of frequencies. It would seem 
reasonable to assume that the dynamics should be 
smoothed less for higher frequencies than for appre- 
ciably lower frequencies. How much attention must 
be paid to this issue? Can one make the crude ap- 
proximation of using smoothing for the center fre- 
quency of a pulse? 
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In pulsed experiments, the early arrivals form 
branches which correspond to wave energy propa- 
gating at the larger angles near 9 max . Depending on 
the range of propagation, the late arrival portion of 
the signal may be constrained to a narrower range of 
propagation angles. Is more smoothing appropriate 
for this portion due to the 6 max factor in A m i„? 

The measure Cav is quite generally a function of 
range. Yet, we mainly used 1000 km propagation in 
our calculations to determine the optimal amount 
of smoothing and neglected the range dependence; 
we did note however a weak range dependence. Re- 
call that several approximations are made arriving 
at the parabolic equation or other one-way, small- 
angle approximation variants. The neglected terms 
may also put range dependence in the propagation, 
and it would not make sense to try to be more accu- 
rate with the smoothing than the level of these other 
approximations. A more detailed understanding of 
the effects of neglected terms would be desirable. 

Although, there is significantly less micro-folding 
for the smoothed than for the unsmoothed poten- 
tials, there is still uncertainty as to how much im- 
provement is gained for the optimal smoothing. This 
could be made clear by carrying out the full de- 
tailed semiclassical construction and comparing it 
to the wave field propagation; we will carry this out 
in Ref. 19. A much deeper understanding would 
come from a full theory based on applying the area- 
(Ao/27r) rule discussed in Sec. IV C. It would give 
the most precise answers possible to questions of 
which structures are physically irrelevant and which 
method removes them in the most optimal way. We 
are pursuing this investigation because only by sep- 
arating out the physically irrelevant fine structures 
can we begin to fully address the ray chaos conun- 
drum and know whether it can be overcome. 
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APPENDIX A: THE SPLIT-OPERATOR, 
FAST FOURIER TRANSFORM METHOD 

The parabolic equation in Eq. (4) describes the 
propagation of an acoustic wave with Hamiltonian 
H = p 2 /2 + V, where p 2 /2 and V denote the kinetic 
and potential energies. A wave field can be advanced 
exactly through the application of the unitary 



propagation operator exp(— iko J Hdr). The split- 
operator Fourier transform method 22 approximates 
this operator using e A+B w e A / 2 e B e A / 2 , where A is 
taken to be -ik J{p 2 /2) dr = -(i/k ) J(k 2 /2) dr 
and B is taken to be —iko J V( z , r)dr. Inserting a 
Fourier transform identity and rearranging terms be- 
fore integrating gives a formula for the propagation 
of a wave field, ^f u (z, r), at a range r to a wave field, 
*&u,(z, r'), at a range r' = r + Ar 



* u (zS) = F- 1 



,A/2p 



e B F~ 



e A ' 2 F[H> u (z,r)] 



(Al) 



where F and F 1 are the forward and backward 
Fourier transforms, respectively. This equation has 

error O ^(fco Ar) 3 ^j due to the operator approxima- 
tion. We approximate the integral V(z, r)dr w 
Ar [V(z, r) + V(z 1 r')}/2 and perform the integration 

J r r '(fc 2 /4) dr = Ar fc 2 /4. 

The wave fields in this paper are calculated over 
a vertical grid of z G [—2,7] km. The reflection 
boundary condition at the surface is not enforced in 
favor of the wave experiencing a soft reflection from 
the potential rather than a hard reflection from the 
surface. Wave energy which is reflected from the 
surface is eventually absorbed by the bottom in long 
range propagations so that this energy is negligible 
at a receiver. The soft reflections of the wave are 
due only to the background portion of the potential 
(Munk's canonical model in Eq. (8)) whose effects 
have been extended above the surface, z < 0. The 
internal wave fluctuations from Eq. (9) are cut off 
by the surface filter in Eq. (10) so that they don't 
have an effect on the wave above the surface. 

The grid size for the propagation is chosen to be 
dependent on the source frequency (to ensure proper 
sampling of the source in the horizontal and verti- 
cal directions) and the maximum number of vertical 
modes, J (to ensure proper sampling of the smallest 
wavelengths in each vertical mode). The grid num- 
ber in the depth direction is purposely chosen to be 
a power of 2 to allow the use of a fast Fourier trans- 
form for the split-operator Fourier method. Specif- 
ically, for the source frequencies 25, 75, 150, 250 Hz, 
we chose Ar = 0.01, 0.01, 0.005, 0.0025 km and Az = 
9/N km where N = 1024,2048,2048,2048, respec- 
tively. These values are large enough to guarantee 
proper convergence of the split-operator method for 
the propagation. 
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APPENDIX B: IMPLEMENTATION OF 
INTERNAL WAVE SOUND SPEED MODEL 

The efficient numerical scheme devised by Colosi 
and Brown 17 generates a random ensemble of in- 
ternal wave effects for the sound speed model, 
Sci w (z,r)/co, through the following equation: 



9 7T 



ScjJ^r) 24.5 2B EAk r 

= #o V -Tj- exp(-3z/2B) (Bl) 



M 



x 



co 

7m.« k max I — 

Yl sin(j<(^))W 2 J , fc " 2 cos(^3, fcr + fc r r) 



j — 1 k r — k m i 

where 



1 

+ - 



1 



V^+i + i 

/? 2 + l ' 2(/32 + 1 )| "M v ^TT-i 



In 



(B2) 



A single random seed generates the random phases, 
4>i,k r € [0j27r), for each internal wave with vertical 
mode, j, and horizontal wavenumber, fc r . These ran- 
dom phases give the ocean a different internal wave 
realization for each random seed. All calculations 
in this paper were done with a single realization of 
the internal wave field, but all results are similar 
for averages over ensemble of internal wave fields as 
well. Each internal wave in the superposition has the 
statistics of the Garrett-Munk spectrum. The full 
Garrett-Munk energy of E = 6.3 x 10 -5 has been 
used in all calculations. Our calculations are done 
for a latitude of 30° so that the inertial frequency 
is fi = 1 cycle per day. The buoyancy profile is 
assumed to have the form N(z) = N a e~ z / B , where 
-/Vq = 1 cycle per 10 min is the buoyancy frequency at 
the surface. We considered the depth of the ocean 
to be H = 5.0 km, even though we extended the 
propagation range to the region [—2,7] km for the 
reasons described in Appendix A. 

The particular functional forms and constants 
used in this paper are as used by Colosi and Brown. 
Some of these forms and constants have already been 
identified in the body of the paper (i.e. near Eq. (9)), 
while the others are listed here. We took the gravita- 
tional acceleration g = 9.81 m/s 2 , M — (717, — l)/2j 2 
and the principle mode number j» = 3. We took 512 
horizontal internal wave numbers equally spaced by 
Ak r for k r G 27r[0.01, 1.0] cycles per km. In the ex- 
pression for Ij.k r , we took kj = fiTrj/N B and the 
ratio (3 — k r /kj. 
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